# create the homophily matrices needed for modeling

#Set Directory
setwd('/Users/drschaef/Google Drive/Projects/Adriana/Raw data/Longitudinal Data/')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Raw data/Longitudinal Data/')

library(network)
library(sna)
library(statnet)
library(dplyr)
library(RSiena)

### M_ data (N=1795 in raw files)
dat <- data.frame(readRDS('attrNoRaceM_220523.rds'),readRDS('attrRaceM_220829_noAMENA.rds')[,-1])
mat1 <-  readRDS("M_W1.network.rds")
mat2 <-  readRDS("M_W2.network.rds")
mat3 <-  readRDS("M_W3.network.rds")
eca1 <- readRDS("ecaM1_220608.rds")
eca2 <- readRDS("ecaM2_220608.rds")
skin <- readRDS("skinM_.rds")

# append skin tone measures to dat
dat <- merge(x=dat, y=skin, all.x=T, all.y=F, by.x='M_D', by.y='MFID')

# missing attribute data kept as NA

##Checking how many having missing WPID for at least 1 wave
table(dat$W1PIDNA)
table(dat$W2PIDNA)
table(dat$W1PIDNA,dat$W2PIDNA)
table(dat$W3PIDNA)

###Indicate which waves of PID is missing
dat$whichPIDmissing <- 0
dat$whichPIDmissing[dat$W1PIDNA==0 & dat$W2PIDNA==0 & dat$W3PIDNA==0] <- "000"
dat$whichPIDmissing[dat$W1PIDNA==1 & dat$W2PIDNA==0 & dat$W3PIDNA==0] <- "100"
dat$whichPIDmissing[dat$W1PIDNA==0 & dat$W2PIDNA==1 & dat$W3PIDNA==0] <- "020"
dat$whichPIDmissing[dat$W1PIDNA==0 & dat$W2PIDNA==0 & dat$W3PIDNA==1] <- "003"
dat$whichPIDmissing[dat$W1PIDNA==1 & dat$W2PIDNA==1 & dat$W3PIDNA==0] <- "120"
dat$whichPIDmissing[dat$W1PIDNA==1 & dat$W2PIDNA==0 & dat$W3PIDNA==1] <- "103"
dat$whichPIDmissing[dat$W1PIDNA==0 & dat$W2PIDNA==1 & dat$W3PIDNA==1] <- "023"
dat$whichPIDmissing[dat$W1PIDNA==1 & dat$W2PIDNA==1 & dat$W3PIDNA==1] <- "123"

#Do some cross tabs for descriptive stats
#table(dat$whichPIDmissing, useNA = "ifany")
#table(dat$Gradecohort, useNA = "ifany")    
#( temp <- table(dat$whichPIDmissing, dat$Gradecohort, useNA = "ifany") )
#table( dat$whichPIDmissing, dat$Gradecohort, (!is.na(dat$chooseOne1) & !is.na(dat$chooseOne2)), useNA = "ifany")

# subset to respondents & valid race at the given period3
# - non-respondents don't have a valid race, and it can't be imputed from other waves
# wave 1 -> 2
keep12 <- (dat$W1PIDNA==0 | dat$W2PIDNA==0) & !is.na(dat$chooseOne1) & !is.na(dat$chooseOne2) & dat$Gradecohort<12
attM12 <- select(filter(dat, keep12,), everything())    # N=667 respondents
matM1_ <- mat1[keep12, keep12]
matM_2 <- mat2[keep12, keep12]

# wave 2 -> 3
keep23 <- (dat$W2PIDNA==0 | dat$W3PIDNA==0) & !is.na(dat$chooseOne2) & !is.na(dat$chooseOne3) & dat$Gradecohort<12
attM23 <- select(filter(dat, keep23,) , everything())    # N=1011 respondents
matM2_ <- mat2[keep23, keep23]
matM_3 <- mat3[keep23, keep23]

# wave 1 -> 3
keep13 <- (dat$W1PIDNA==0 | dat$W3PIDNA==0) & !is.na(dat$chooseOne1) & !is.na(dat$chooseOne3) & dat$Gradecohort<12
attM13 <- select(filter(dat, keep13,) , everything())    # N=616 respondents
matM13_ <- mat1[keep13, keep13]
matM_13 <- mat3[keep13, keep13]

# subset ECAs and create overlap matrices
# add fake second matrix (used when keeping ECA exogenous but testing effects on race )
temp <- as.matrix(eca1[keep12,-1])
temp[is.na(temp)|temp==10] <- 0
ecaM12 <- temp %*% t(temp)
ecaM12dum <- ecaM12
ecaM12dum[ecaM12 > 0] <- 1
	ecaRigged <- ecaM12dum
	second0 <- which(ecaRigged==0)[2]
	second1 <- which(ecaRigged==1)[2]
	ecaRigged[second0] <- 1
	ecaRigged[second1] <- 0
ecaM12list <- list(ecaM12dum, ecaRigged)

temp <- as.matrix(eca2[keep23,-1])
temp[is.na(temp)|temp==10] <- 0
ecaM23 <- temp %*% t(temp)
ecaM23dum <- ecaM23
ecaM23dum[ecaM23 > 0] <- 1
	ecaRigged <- ecaM23dum
	second0 <- which(ecaRigged==0)[2]
	second1 <- which(ecaRigged==1)[2]
	ecaRigged[second0] <- 1
	ecaRigged[second1] <- 0
ecaM23list <- list(ecaM23dum, ecaRigged)

# for time 1 to time 3 use the real data
temp <- as.matrix(eca1[keep13,-1])  # wave 1
temp[is.na(temp)|temp==10] <- 0
ecaM13_ <- temp %*% t(temp)
ecaM13_dum <- ecaM13_
ecaM13_dum[ecaM13_ > 0] <- 1
temp <- as.matrix(eca2[keep13,-1])  # wave 3
temp[is.na(temp)|temp==10] <- 0
ecaM_13 <- temp %*% t(temp)
ecaM_13dum <- ecaM_13
ecaM_13dum[ecaM_13 > 0] <- 1
ecaM13list <- list(ecaM13_dum, ecaM_13dum)


# create two mode current race matrices
# wave 1 -> 2
# - start by creating structural zeros based on heritage
raceM1_ <- as.matrix(attM12[,which(names(attM12)=='A'): which(names(attM12)=='MULTI')], ncol=7)
raceM1_[raceM1_==0] <- 10
raceM1_[raceM1_==1] <- 0
raceM1_[attM12$multi==1,7] <- 0   # anyone with multi heritage could choose multi
raceM_2 <- raceM1_
#  - plug in current race for each wave
for (i in 1:dim(attM12)[1]) {
	if (attM12$chooseOne1[i] <= 7) { raceM1_[i, attM12$chooseOne1[i]] <- (1+attM12$mono[i]*10) }  # monos become 11, multis become 1
	else if (attM12$chooseOne1[i] > 7) { raceM1_[i, (attM12$chooseOne1[i]-2)] <- (1+attM12$mono[i]*10) }  # other & multi shift down
	if (attM12$chooseOne2[i] <= 7) { raceM_2[i, attM12$chooseOne2[i]] <- (1+attM12$mono[i]*10) }
	else if (attM12$chooseOne2[i] > 7) { raceM_2[i, (attM12$chooseOne2[i]-2)] <- (1+attM12$mono[i]*10) }   # other & multi shift down
	}

# wave 2 -> 3
# - start by creating structural zeros based on heritage
raceM2_ <- as.matrix(attM23[,which(names(attM23)=='A'): which(names(attM23)=='MULTI')], ncol=7)
raceM2_[raceM2_==0] <- 10
raceM2_[raceM2_==1] <- 0
raceM2_[attM23$multi==1,7] <- 0   # anyone with multi heritage could choose multi
raceM_3 <- raceM2_
#  - plug in current race for each wave
for (i in 1:dim(attM23)[1]) {
	if (attM23$chooseOne2[i] <= 7) { raceM2_[i, attM23$chooseOne2[i]] <- (1+attM23$mono[i]*10) }
	else if (attM23$chooseOne2[i] > 7) { raceM2_[i, (attM23$chooseOne2[i]-2)] <- (1+attM23$mono[i]*10) }  # other & multi shift down
	if (attM23$chooseOne3[i] <= 7) { raceM_3[i, attM23$chooseOne3[i]] <- (1+attM23$mono[i]*10) }
	else if (attM23$chooseOne3[i] > 7) { raceM_3[i, (attM23$chooseOne3[i]-2)] <- (1+attM23$mono[i]*10) }   # other & multi shift down
	}

# wave 1 -> 3
# - start by creating structural zeros based on heritage
raceM13_ <- as.matrix(attM13[,which(names(attM13)=='A'): which(names(attM13)=='MULTI')], ncol=7)
raceM13_[raceM13_==0] <- 10
raceM13_[raceM13_==1] <- 0
raceM13_[attM13$multi==1,7] <- 0   # anyone with multi heritage could choose multi
raceM_13 <- raceM13_
#  - plug in current race for each wave
for (i in 1:dim(attM13)[1]) {
	if (attM13$chooseOne1[i] <= 7) { raceM13_[i, attM13$chooseOne1[i]] <- (1+attM13$mono[i]*10) }
	else if (attM13$chooseOne1[i] > 7) { raceM13_[i, (attM13$chooseOne1[i]-2)] <- (1+attM13$mono[i]*10) }  # other & multi shift down
	if (attM13$chooseOne3[i] <= 7) { raceM_13[i, attM13$chooseOne3[i]] <- (1+attM13$mono[i]*10) }
	else if (attM13$chooseOne3[i] > 7) { raceM_13[i, (attM13$chooseOne3[i]-2)] <- (1+attM13$mono[i]*10) }   # other & multi shift down
	}
	
r.c.m12 <- table(attM12$chooseOne1,attM12$chooseOne2)
r.c.m23 <- table(attM12$chooseOne2,attM12$chooseOne3)

### S_ data (N=3298 in raw files) 
dat <- data.frame(readRDS('attrNoRaceS_220523.rds'),readRDS('attrRaceS_220829_noAMENA.rds')[,-1])
mat1 <-  readRDS("S_W1.network.rds")
mat2 <-  readRDS("S_W2.network.rds")
mat3 <-  readRDS("S_W3.network.rds")
eca1 <- readRDS("ecaS1_220608.rds")
eca2 <- readRDS("ecaS2_220608.rds")
skin <- readRDS("skinS_.rds")

# append skin tone measures to dat
dat <- merge(x=dat, y=skin, all.x=T, all.y=F, by.x='M_D', by.y='MFID')

# missing attribute data kept as NA

##Checking how many having missing WPID for at least 1 wave
table(dat$W1PIDNA)
table(dat$W2PIDNA)
table(dat$W3PIDNA)

###Indicate which waves of PID is missing
dat$whichPIDmissing <- 0
dat$whichPIDmissing[dat$W1PIDNA==0 & dat$W2PIDNA==0 & dat$W3PIDNA==0] <- "000"
dat$whichPIDmissing[dat$W1PIDNA==1 & dat$W2PIDNA==0 & dat$W3PIDNA==0] <- "100"
dat$whichPIDmissing[dat$W1PIDNA==0 & dat$W2PIDNA==1 & dat$W3PIDNA==0] <- "020"
dat$whichPIDmissing[dat$W1PIDNA==0 & dat$W2PIDNA==0 & dat$W3PIDNA==1] <- "003"
dat$whichPIDmissing[dat$W1PIDNA==1 & dat$W2PIDNA==1 & dat$W3PIDNA==0] <- "120"
dat$whichPIDmissing[dat$W1PIDNA==1 & dat$W2PIDNA==0 & dat$W3PIDNA==1] <- "103"
dat$whichPIDmissing[dat$W1PIDNA==0 & dat$W2PIDNA==1 & dat$W3PIDNA==1] <- "023"
dat$whichPIDmissing[dat$W1PIDNA==1 & dat$W2PIDNA==1 & dat$W3PIDNA==1] <- "123"

#Do some cross tabs for descriptive stats
table(dat$whichPIDmissing, useNA = "ifany")
table(dat$Gradecohort, useNA = "ifany")    # 52 missing, but they are dropped below
table(dat$whichPIDmissing, dat$Gradecohort, useNA = "ifany")

# subset to respondents & valid race at the given period3
# - non-respondents don't have a valid race, and it can't be imputed from other waves
# wave 1 -> 2
keep12 <- (dat$W1PIDNA==0 | dat$W2PIDNA==0) & !is.na(dat$chooseOne1) & !is.na(dat$chooseOne2) & dat$Gradecohort<12
attS12 <- select(filter(dat, keep12,), everything())    # N=1251 respondents
matS1_ <- mat1[keep12, keep12]
matS_2 <- mat2[keep12, keep12]

# wave 2 -> 3
keep23 <- (dat$W2PIDNA==0 | dat$W3PIDNA==0) & !is.na(dat$chooseOne2) & !is.na(dat$chooseOne3) & dat$Gradecohort<12
attS23 <- select(filter(dat, keep23,) , everything())    # N=1559 respondents
matS2_ <- mat2[keep23, keep23]
matS_3 <- mat3[keep23, keep23]

# wave 1 -> 3
keep13 <- (dat$W1PIDNA==0 | dat$W3PIDNA==0) & !is.na(dat$chooseOne1) & !is.na(dat$chooseOne3) & dat$Gradecohort<12
attS13 <- select(filter(dat, keep13,) , everything())    # N=1010 respondents
matS13_ <- mat1[keep13, keep13]
matS_13 <- mat3[keep13, keep13]

# subset ECAs and create overlap matrices
# add fake second matrix (used when keeping ECA exogenous but testing effects on race )
temp <- as.matrix(eca1[keep12,-1])
temp[is.na(temp)|temp==10] <- 0
ecaS12 <- temp %*% t(temp)
ecaS12dum <- ecaS12
ecaS12dum[ecaS12 > 0] <- 1
	ecaRigged <- ecaS12dum
	second0 <- which(ecaRigged==0)[2]
	second1 <- which(ecaRigged==1)[2]
	ecaRigged[second0] <- 1
	ecaRigged[second1] <- 0
ecaS12list <- list(ecaS12dum, ecaRigged)

temp <- as.matrix(eca2[keep23,-1])
temp[is.na(temp)|temp==10] <- 0
ecaS23 <- temp %*% t(temp)
ecaS23dum <- ecaS23
ecaS23dum[ecaS23 > 0] <- 1
	ecaRigged <- ecaS23dum
	second0 <- which(ecaRigged==0)[2]
	second1 <- which(ecaRigged==1)[2]
	ecaRigged[second0] <- 1
	ecaRigged[second1] <- 0
ecaS23list <- list(ecaS23dum, ecaRigged)

# for time 1 to time 3 use the real data
temp <- as.matrix(eca1[keep13,-1])  # wave 1
temp[is.na(temp)|temp==10] <- 0
ecaS13_ <- temp %*% t(temp)
ecaS13_dum <- ecaS13_
ecaS13_dum[ecaS13_ > 0] <- 1
temp <- as.matrix(eca2[keep13,-1])  # wave 3
temp[is.na(temp)|temp==10] <- 0
ecaS_13 <- temp %*% t(temp)
ecaS_13dum <- ecaS_13
ecaS_13dum[ecaS_13 > 0] <- 1
ecaS13list <- list(ecaS13_dum, ecaS_13dum)


# create two mode current race matrices
# wave 1 -> 2
# - start by creating structural zeros based on heritage
raceS1_ <- as.matrix(attS12[,which(names(attS12)=='A'): which(names(attS12)=='MULTI')], ncol=7)
raceS1_[raceS1_==0] <- 10
raceS1_[raceS1_==1] <- 0
raceS1_[attS12$multi==1,7] <- 0   # anyone with multi heritage could choose multi
raceS_2 <- raceS1_
#  - plug in current race for each wave
for (i in 1:dim(attS12)[1]) {
	if (attS12$chooseOne1[i] <= 7) { raceS1_[i, attS12$chooseOne1[i]] <- (1+attS12$mono[i]*10) }
	else if (attS12$chooseOne1[i] > 7) { raceS1_[i, (attS12$chooseOne1[i]-2)] <- (1+attS12$mono[i]*10) }  # other & multi shift down
	if (attS12$chooseOne2[i] <= 7) { raceS_2[i, attS12$chooseOne2[i]] <- (1+attS12$mono[i]*10) }
	else if (attS12$chooseOne2[i] > 7) { raceS_2[i, (attS12$chooseOne2[i]-2)] <- (1+attS12$mono[i]*10) }   # other & multi shift down
	}

# wave 2 -> 3
# - start by creating structural zeros based on heritage
raceS2_ <- as.matrix(attS23[,which(names(attS23)=='A'): which(names(attS23)=='MULTI')], ncol=7)
raceS2_[raceS2_==0] <- 10
raceS2_[raceS2_==1] <- 0
raceS2_[attM23$multi==1,7] <- 0   # anyone with multi heritage could choose multi
raceS_3 <- raceS2_
#  - plug in current race for each wave
for (i in 1:dim(attS23)[1]) {
	if (attS23$chooseOne2[i] <= 7) { raceS2_[i, attS23$chooseOne2[i]] <- (1+attS23$mono[i]*10) }
	else if (attS23$chooseOne2[i] > 7) { raceS2_[i, (attS23$chooseOne2[i]-2)] <- (1+attS23$mono[i]*10) }  # other & multi shift down
	if (attS23$chooseOne3[i] <= 7) { raceS_3[i, attS23$chooseOne3[i]] <- (1+attS23$mono[i]*10) }
	else if (attS23$chooseOne3[i] > 7) { raceS_3[i, (attS23$chooseOne3[i]-2)] <- (1+attS23$mono[i]*10) }   # other & multi shift down
	}

# wave 1 -> 3
# - start by creating structural zeros based on heritage
raceS13_ <- as.matrix(attS13[,which(names(attS13)=='A'): which(names(attS13)=='MULTI')], ncol=7)
raceS13_[raceS13_==0] <- 10
raceS13_[raceS13_==1] <- 0
raceS13_[attM13$multi==1,7] <- 0   # anyone with multi heritage could choose multi
raceS_13 <- raceS13_
#  - plug in current race for each wave
for (i in 1:dim(attS13)[1]) {
	if (attS13$chooseOne1[i] <= 7) { raceS13_[i, attS13$chooseOne1[i]] <- (1+attS13$mono[i]*10) }
	else if (attS13$chooseOne1[i] > 7) { raceS13_[i, (attS13$chooseOne1[i]-2)] <- (1+attS13$mono[i]*10) }  # other & multi shift down
	if (attS13$chooseOne3[i] <= 7) { raceS_13[i, attS13$chooseOne3[i]] <- (1+attS13$mono[i]*10) }
	else if (attS13$chooseOne3[i] > 7) { raceS_13[i, (attS13$chooseOne3[i]-2)] <- (1+attS13$mono[i]*10) }   # other & multi shift down
	}

r.c.s12 <- table(attS12$chooseOne1,attS12$chooseOne2)
r.c.s23 <- table(attS12$chooseOne2,attS12$chooseOne3)
# race change: 1=A, 2=B, 3=L, 4=W, 5=N, 8=O, 9=multi
r.c.m12 + r.c.m23  # M_ 
r.c.s12 + r.c.s23  # S_ 
(temp <- r.c.m23 + r.c.s23 + r.c.m12  + r.c.s12)  # full sample
sum(temp) - sum(diag(temp))  #258 change -> became 186 on 10/14??

# size of each ethnoracial group in schools (ever claimed identity)
pPopM <- c(.297,.25,.114,.52,.064,.061,.043)  # multi could be recoded as .241 (% who could claim multi)
pPopS <- c(.09,.337,.395,.501,.156,.023,.032) # multi could be recoded as .379 (% who could claim multi)

#save.image("~/Google Drive/Projects/Adriana/Papers/Homophily/allDatPreSIENA221014.RData")
load("~/Google Drive/Projects/Adriana/Papers/Homophily/allDatPreSIENA221014.RData")

# function to create heritage homophily matrices
makeHeritageSim <- function(ATT) {
	rDums <- as.matrix(ATT[,which(names(ATT)=='A'): which(names(ATT)=='N')], ncol=5)
	nInCommon <- rDums %*% t(rDums)
	N <- dim(ATT)[1]
	nRace <- rowSums(rDums)
	maxInCommon <- pmax(matrix(rep(nRace,N),ncol=N), t(matrix(rep(nRace,N),ncol=N)) )
	maxInCommon[maxInCommon==0] <- .000001
	homophFull <- nInCommon / maxInCommon
	homophPart <- homophFull
	homophPart[homophPart==1] <- 0
	homophPart[homophPart>0] <- 1
	homophFull[homophFull<1] <- 0
	diag(homophFull) <- 0
	diag(homophPart) <- 0
	# separate into multi-multi dyads and multi-mono dyads
	momoFull <- homophFull * outer(ATT$mono==1, ATT$mono==1, '*')
	mumuFull <- homophFull * outer(ATT$multi==1, ATT$multi==1, '*')
	mumuPart <- homophPart * outer(ATT$multi==1, ATT$multi==1, '*')
	mumoPart <- homophPart * (outer(ATT$mono==1, ATT$multi==1, '*') + 
							 outer(ATT$multi==1, ATT$mono==1, '*'))
	mumu <- (homophFull + homophPart) * outer(ATT$multi==1, ATT$multi==1, '*')
	mumo <- homophPart * outer(ATT$multi==1, ATT$mono==1, '*')
	momu <- homophPart * outer(ATT$mono==1, ATT$multi==1, '*')
	return(array(c(homophFull, homophPart, momoFull, mumuFull, mumuPart, mumoPart, mumu, mumo, momu), dim=c(N,N,9)))
	}

# a function to create the SIENA object
# as part of the function, create heritage homophily dyadic covariate (scaled using .6)
#  - 	HOM <- HOM1[,,1] + .6*HOM1[,,2] # uses hSim layers 1 and 2

makeSienaData <- function (NET1, NET2, ATT, ECA, RACE1, RACE2, POPSIZE, STARTWAVE, DVS=NULL) {
	# calculate homophily on heritage
	hSim <- makeHeritageSim(ATT)
	erSimScaled <- hSim[,,1] + .6*hSim[,,2]
	erSimNomomo <- hSim[,,4] + hSim[,,5] + hSim[,,6]

	# counts
	nActors <- dim(NET1)[1]
	nRace <- dim(RACE1)[2]
	
	# define different node sets:
	students <- sienaNodeSet(nActors, nodeSetName="students")
	erGroups <- sienaNodeSet(nRace, nodeSetName="erGroups")

	# dependent networks
	friendship <- sienaNet(array( c(NET1, NET2), dim = c(nActors, nActors, 2)), nodeSet="students")
	raceNet <- sienaNet(array( c(RACE1, RACE2), dim = c(nActors, nRace, 2)), "bipartite", nodeSet=c("students","erGroups"))

	# endogenous version of ECA (rigged to be exogenous)
	ecaEndog <- sienaNet(array( c(ECA[[1]], ECA[[2]]), dim = c(nActors, nActors, 2)), nodeSet="students")

	# student attributes
	female <- coCovar(ATT$Female, nodeSet="students")
	gradecohort <- coCovar(ATT$Gradecohort, nodeSet="students")
	pared <- coCovar(ATT$ParEd, nodeSet="students")
	pared8 <- coCovar(ATT$pared8, nodeSet="students")
	immgen <- coCovar(ATT$ImmigGen, nodeSet="students")
	multi <- coCovar(ATT$multi, center=F, nodeSet="students")
	mono <- coCovar(ATT$mono, center=F, nodeSet="students")
	nis <- coCovar(ATT$NIS, nodeSet="students")
	perla <- coCovar(ATT$PERLA, nodeSet="students")
	racedif <- coCovar(ATT$RACEDIF, nodeSet="students")
	egoA <- coCovar(ATT$A, nodeSet="students")
	egoB <- coCovar(ATT$B, nodeSet="students")
	egoL <- coCovar(ATT$L, nodeSet="students")
	egoW <- coCovar(ATT$W, nodeSet="students")
	egoN <- coCovar(ATT$N, nodeSet="students")
	egoO <- coCovar(ATT$O, nodeSet="students")
	egoM <- coCovar(ATT$MUL, nodeSet="students")
	
	# exogenously changing covariates (treated as constant within multigroup setup)
	# - academic performance
	tempAcad <- ATT[,grepl("academic", names(ATT), fixed=T)]
	tempAcad[tempAcad==-9] <- NA
	tempAcad[is.na(tempAcad[,1]),1] <- tempAcad[is.na(tempAcad[,1]),2]  # if t1 missing, impute with t2
	acad <- coCovar(tempAcad[,STARTWAVE], nodeSet="students")
	tempAcad9 <- ATT[,grepl("acad9", names(ATT), fixed=T)]
	tempAcad9[is.na(tempAcad9[,1]),1] <- tempAcad9[is.na(tempAcad9[,1]),2]  # if t1 missing, impute with t2
	acad9 <- coCovar(tempAcad9[,STARTWAVE], nodeSet="students")
	# - exogeneous version of current race
	tempRace <- data.frame(ATT$chooseOne1, ATT$chooseOne2, ATT$chooseOne3)
	raceExog <- coCovar(tempRace[,STARTWAVE], nodeSet="students")
		
	# - discrimination, overall
	temp5 <- data.frame(ATT$disc1, ATT$disc2, ATT$disc3)
	disc <- coCovar(temp5[,STARTWAVE], nodeSet="students")
	
	# race attributes: other = 6 (reference category)
	asian <- coCovar(c(1,0,0,0,0,0,0), center=F, nodeSet='erGroups')
	black <- coCovar(c(0,1,0,0,0,0,0), center=F, nodeSet='erGroups')
	latin <- coCovar(c(0,0,1,0,0,0,0), center=F, nodeSet='erGroups')
	white <- coCovar(c(0,0,0,1,0,0,0), center=F, nodeSet='erGroups')
	natam <- coCovar(c(0,0,0,0,1,0,0), center=F, nodeSet='erGroups')
	orace <- coCovar(c(0,0,0,0,0,1,0), center=F, nodeSet='erGroups')
	mrace <- coCovar(c(0,0,0,0,0,0,1), center=F, nodeSet='erGroups')
	pPop <- coCovar(POPSIZE, nodeSet='erGroups')

	# dyadic covariates (student X student)
	ecaDum <- coDyadCovar(ECA[[1]], center=F, nodeSets=c("students", "students"))
	herRaceSim <- coDyadCovar(erSimScaled, center=F, nodeSets=c("students", "students"))
	herRaceSim2 <- coDyadCovar(erSimNomomo, center=F, nodeSets=c("students", "students"))
	mumu <- coDyadCovar(hSim[,,7], center=F, nodeSets=c("students", "students"))
	mumo <- coDyadCovar(hSim[,,8], center=F, nodeSets=c("students", "students"))
	momu <- coDyadCovar(hSim[,,9], center=F, nodeSets=c("students", "students"))

	# dyadic covariates (student X race)
	immgenA <- coDyadCovar(outer(immgen, asian, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenB <- coDyadCovar(outer(immgen, black, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenL <- coDyadCovar(outer(immgen, latin, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenW <- coDyadCovar(outer(immgen, white, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenN <- coDyadCovar(outer(immgen, natam, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenO <- coDyadCovar(outer(immgen, orace, "*"), center=T, nodeSets=c("students","erGroups"))
	immgenM <- coDyadCovar(outer(immgen, mrace, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaA <- coDyadCovar(outer(perla, asian, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaB <- coDyadCovar(outer(perla, black, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaL <- coDyadCovar(outer(perla, latin, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaW <- coDyadCovar(outer(perla, white, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaN <- coDyadCovar(outer(perla, natam, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaO <- coDyadCovar(outer(perla, orace, "*"), center=T, nodeSets=c("students","erGroups"))
	perlaM <- coDyadCovar(outer(perla, mrace, "*"), center=T, nodeSets=c("students","erGroups"))
	discA <- coDyadCovar(outer(temp5[,STARTWAVE], asian, "*"), center=T, nodeSets=c("students","erGroups"))
	discB <- coDyadCovar(outer(temp5[,STARTWAVE], black, "*"), center=T, nodeSets=c("students","erGroups"))
	discL <- coDyadCovar(outer(temp5[,STARTWAVE], latin, "*"), center=T, nodeSets=c("students","erGroups"))
	discW <- coDyadCovar(outer(temp5[,STARTWAVE], white, "*"), center=T, nodeSets=c("students","erGroups"))
	discN <- coDyadCovar(outer(temp5[,STARTWAVE], natam, "*"), center=T, nodeSets=c("students","erGroups"))
	discO <- coDyadCovar(outer(temp5[,STARTWAVE], orace, "*"), center=T, nodeSets=c("students","erGroups"))
	discM <- coDyadCovar(outer(temp5[,STARTWAVE], mrace, "*"), center=T, nodeSets=c("students","erGroups"))


	myRaceFriendsData <- sienaDataCreate( friendship, raceNet,
	       female, gradecohort, pared, pared8, immgen, acad, acad9, multi, mono, nis, perla, racedif,
	       asian, black, latin, white, natam, orace, mrace, pPop,
	       disc, egoA, egoB, egoL, egoW, egoN, egoO, egoM,
	       ecaDum, herRaceSim, herRaceSim2, mumu, mumo, momu, 
	       nodeSets=list(students, erGroups))
	
	myRaceFriendsDataRI <- sienaDataCreate( friendship, raceNet,
	       female, gradecohort, pared, pared8, immgen, acad, acad9, multi, mono, nis, perla, racedif,
	       asian, black, latin, white, natam, orace, mrace, pPop,
	       disc, egoA, egoB, egoL, egoW, egoN, egoO, egoM,
	       ecaDum, herRaceSim, herRaceSim2, mumu, mumo, momu,  
				  immgenA, immgenB, immgenL, immgenN, immgenW, immgenO, immgenM,
				  perlaA, perlaB, perlaL, perlaN, perlaW, perlaO, perlaM,
				  discA, discB, discL, discN, discW, discO, discM,
	       nodeSets=list(students, erGroups))
	       
	RET <- myRaceEvolutionData 	
	if (is.null(DVS)) { RET <- myCoEvolutionData }
	else if (DVS=='race+friendship') { RET <- myRaceFriendsData }
	else if (DVS=='race+fr_ri') { RET <- myRaceFriendsDataRI }
		return(RET)
	
	}


# race + friendship
rfnetM12 <- makeSienaData(matM1_, matM_2, attM12, ecaM12list, raceM1_, raceM_2, pPopM, 1, 'race+friendship')
rfnetM23 <- makeSienaData(matM2_, matM_3, attM23, ecaM23list, raceM2_, raceM_3, pPopM, 2, 'race+friendship')
rfnetS12 <- makeSienaData(matS1_, matS_2, attS12, ecaS12list, raceS1_, raceS_2, pPopS, 1, 'race+friendship')
rfnetS23 <- makeSienaData(matS2_, matS_3, attS23, ecaS23list, raceS2_, raceS_3, pPopS, 2, 'race+friendship')
rfnetM123 <- sienaGroupCreate(list(rfnetM12, rfnetM23))
rfnetS123 <- sienaGroupCreate(list(rfnetS12, rfnetS23))

# race + friendship for Relative Importance analysis
rfnetM12ri <- makeSienaData(matM1_, matM_2, attM12, ecaM12list, raceM1_, raceM_2, pPopM, 1, 'race+fr_ri')
rfnetM23ri <- makeSienaData(matM2_, matM_3, attM23, ecaM23list, raceM2_, raceM_3, pPopM, 2, 'race+fr_ri')
rfnetS12ri <- makeSienaData(matS1_, matS_2, attS12, ecaS12list, raceS1_, raceS_2, pPopS, 1, 'race+fr_ri')
rfnetS23ri <- makeSienaData(matS2_, matS_3, attS23, ecaS23list, raceS2_, raceS_3, pPopS, 2, 'race+fr_ri')
rfnetM123ri <- sienaGroupCreate(list(rfnetM12ri, rfnetM23ri))
rfnetS123ri <- sienaGroupCreate(list(rfnetS12ri, rfnetS23ri))


setwd('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Papers/Homophily/Output')

saveRDS(rfnetM123, file = 'rfnetM123_221014.rds')
saveRDS(rfnetS123, file = 'rfnetS123_221014.rds')
saveRDS(rfnetM123ri, file = 'rfnetM123_221014ri.rds')
saveRDS(rfnetS123ri, file = 'rfnetS123_221014ri.rds')
saveRDS(rfnetM12ri, file = 'rfnetM12_221014ri.rds')
saveRDS(rfnetM23ri, file = 'rfnetM23_221014ri.rds')
saveRDS(rfnetS12ri, file = 'rfnetS12_221014ri.rds')
saveRDS(rfnetS23ri, file = 'rfnetS23_221014ri.rds')


